function PluvioScalogram(sigTot,prDays,neDays,nPerDay,d1,d2,sizeIm,comPart,chPlot)
% PluvioScalogram - Scalogram computation for a rainfall time series
%
% SigSegmData = PluvioScalogram(tSeries,prDays,nrDays,nPerDay,d1,d2,sizeIm,comPart,chPlot)
%
% Let tSeries = [t Pluvio] be a time series of pluviometric data, where
% Pluvio(k) is the rainfall value (in mm) corresponding to the time t(k) 
% espressed in MATLAB serial number form, with t(k+1) - t(k) = 1/24, i.e. 
% 1h, for eack k. Therefore, the sampling frequency is Fs = 24.
% Another requirement about the time series is that t(1) and t(end) are 
% integer numbers, i.e. are dates related to 0h.
% The elements of the input Pluvio can be differential or cumulative values 
% (at some days or fraction of day), according to the user's preference. 
% This function computes the scalogram by means of Continuous Wavelet 
% Transform (CWT). The parameters for CWT computation are:
%   - analytic Morlet wavelet;
%   - 24 voices per octave. 
% For computation efficiency purposes, the CWT is computed by means of 
% preliminary CWT filter bank generation.  
% Each scalogram is computed for the reference date t(k) on the basis of 
% these input data:
%
%   - prDays    : the starting data of the scalogram related to t(k) is
%                 t(k)-prDays, where prDays is expressed in days, under 
%                 the condition that t(k)-prDays >= t(1) (if this condition 
%                 is not fullfilled, the corresponding scalogram is not 
%                 computed). 
%                 If prDays is undefined or empty, the default choice
%                 prDays = 15 is used.
%
%   - neDays    : the ending data of the scalogram related to t(k) is
%                 t(k)+neDays, where nextDays is expressed in days, under 
%                 the condition that t(k)+neDays <= t(end) (if such a
%                 condition is not fullfilled, the corresponding scalogram 
%                 is not computed). 
%                 If neDays is undefined or empty, the default choice
%                 neDays = 0 is used.
%
%   - nPerDay   : number of scalograms per day. The allowed values are
%                 1 (one scalogram per day at 0h), 2 (two scalograms at 
%                 0h and 12 h), 3 (0h, 8h, 16h), 4 (0h, 6h, 12h, 18h, 6
%                 (0h, 4h, 8h, 12h, 16h, 18h, 20h) 12 (even hours), 24.    
%                 If nPerDay is undefined, empty or not an allowed value, 
%                 nPerDay = 1 is used.  
%
%   - d1        : initial date for which the scalograms are computed (this
%                 date must be expressed as serial number). The start time 
%                 of the first computed scalogram is max(d1-prevDays,t(1)). 
%                 If d1 is undefined of empty, d1 = t(1).
%
%   - d2        : final date for which the scalograms are computed (this
%                 date must be expressed as serial number). The end time 
%                 of the last computed scalogram is min(d2+nextDays,t(end)). 
%                 If d2 is undefined of empty, d2 = t(end).
%
%   - sizeIm    : size of the output images. 
%                 If sizeIm is a vector, it is directly used (the number of
%                 channels is unnecessary because the fact that RGB images 
%                 are generated is understood. 
%                 If sizeIm is a scalar, sizeIm = [sizeIm sizeIm] is used.
%                 If sizeIm is undefined or empty, sizeIm = [224 224] is 
%                 used (such a size is the standard one for VGG16, VGG19,
%                 GoogLeNet and ResNet convolutional neural networks). 
%
%   - comPart:    common part of the output filename (folder name, constant 
%                 part of the filename). Each filename is completed with the 
%                 date string. For example, if comPart is the string 
%                       'pluvioData\LongBeach15days', 
%                 the filename of the scalogram of data related to the 15 
%                 days until the date 15 march 2020, 12:00:00 is the string
%                    'pluvioData\LongBeach15days_15-Mar-2020-120000.jpg'.
%                 If comPart is undefined or empty, comPart = '' is used.
%                 If comPart contains a folder path and this folder does
%                 not exist, it is generated.
%
%   - chPlot:     vector of dates for which, besides the file generation, 
%                 the scalogram is also shown in a figure. The dates in
%                 chPlot must be coherent with the choice about nPerDay.
%                 If chPlot is undefined or empty, no figures are shown.
%
% The duration of each signal segment is prDays+neDays. The corresponding 
% length is (prDays+neDays)*Fs+1.
%
% The output SigSegmData is a 4-column cell array such that SigSegmData{k,1} 
% is is the reference time of the k-th row, SigSegmData{k,2} is the initial
% time of the rainfall time series, SigSegmData{k,3} is the corresponding
% end time, SigSegmData{k,4} is the generated filename. In this way, the data
% about the signal segment with contributes to a figure can be stored.
%
% See also PluVelScalogram, trendClass.
%
% SigSegmData = PluvioScalogram(tSeries,prDays,neDays,nPerDay,d1,d2,sizeIm,comPart,chPlot)

% G. Teza, 2020

if nargin < 9
    chPlot = [];
end

if nargin < 8 || isempty(comPart)
    comPart = '';
end
[pathFold,~,~] = fileparts(comPart);
exist(pathFold,'dir')
if exist(pathFold,'dir') ~= 7
    mkdir(pathFold)
end

if nargin < 7 || isempty(sizeIm)
    sizeIm = [224 224];
end
if isscalar(sizeIm)
    sizeIm = [sizeIm sizeIm];
end

if nargin < 6
    d2 = [];
end

if nargin < 5
    d1 = [];
end

if (nargin < 4)||isempty(nPerDay)||~ismember(nPerDay,[1 2 3 4 6 12 24])
    nPerDay = 1;
end

if (nargin < 3)||isempty(neDays)
    neDays = 0;
end

if (nargin < 2)||isempty(prDays)
    prDays = 15;
end

t = sigTot(:,1);
dt = median(diff(t));
if abs(dt-1/24) > 10^6*eps
    error('Invalid input data: the sampling must be 1/24 h');
end
Fs = 24;
sig = sigTot(:,2);
lSigSegm = round((prDays+neDays)*Fs)+1;     % length of each signal segment

if isempty(d1) || (d1 < t(1))
    n1 = 1;
else
    d1 = floor(d1);
    [~,n1] = min(abs(t-d1));
end
if isempty(d2) || (d2 > t(end))
    n2 = numel(t);
else
    d2 = ceil(d2);
    [~,n2] = min(abs(t-d2));
end

nStep = round(24/nPerDay);

% CWT filter bank;
cwtfb = cwtfilterbank(...
    'SignalLength',lSigSegm,...
    'Wavelet','amor',...
    'VoicesPerOctave',24,...
    'samplingFrequency',Fs);

nt = n2-n1+1;

hwb = waitbar(0,'Scalogram generations in progress...');
cmap = jet(128);
nplot = 1;
for k = n1:nStep:n2
    trefk = t(k);
    t1k = trefk-prDays;
    t2k = trefk+neDays; 
    if (t1k >= t(1))&&(t2k <= t(end))
        [~,n1k] = min(abs(t-t1k));
        [~,n2k] = min(abs(t-t2k));
        tk = t(n1k:n2k);
        sigk = sig(n1k:n2k);
        if nonanTS(sigk) && (numel(sigk) == lSigSegm)  % last check
            % Computation by using the precomputed filter bank:
            [cfsk,fk] = cwtfb.wt(sigk);
            % Plot (if required):
            if ~isempty(chPlot) && (nplot <= numel(chPlot))
                if abs(t2k-chPlot(nplot)) < 10^3*eps
                    hfk = figure;
                    pcolor(tk,log2(fk),abs(cfsk));
                    hfk.Colormap = cmap;
                    cl = colorbar;
                    cl.Label.String = 'magnitude';
                    axis tight
                    title(sprintf('Scalogram - from %s to %s',...
                        datestr(t1k),datestr(t2k)));
                    xlabel('\ittime \rm(days)');
                    ylabel('log_2(\itfrequency\rm(cpd))');
                    nplot = nplot+1;
                end
            end
            % Image generation ad saving:
            filenak = [comPart '_' datestr(t2k) '.jpg'];
            filenak = strrep(filenak,' ','-');  % to change " " to "-"  
            filenak = strrep(filenak,':','');   % to remove ":" 
            imk = ind2rgb(im2uint8(rescale(abs(cfsk))),cmap);
            imk = imresize(imk,sizeIm);
            imwrite(imk,filenak);
            waitbar((k-n1)/nt,hwb);
        end
    end
end
close(hwb);

%% Ancillary functions:

function Ina = nonanTS(TS)
% Ina is 1 if TS doeas not contain any NaN or Inf value, 0 elsewhere. 
I1 = isnan(TS);
I2 = isinf(TS);
if (sum(I1) > 0) || (sum(I2) > 0)
    Ina = 0;
else
    Ina = 1;
end
Ina = logical(Ina);